stLearn ：空间轨迹推断

作者 |  周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树。

生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。

stLearn: integrating spatial location, tissue morphology and gene expression to find cell types, cell-cell interactions and spatial trajectories within undissociated tissues

空间转录组学是一种新兴的技术，它将空间信息和组织形态以及表达量融合成空间表达数据。整合这三种类型的数据极大地丰富了人们的想象力，寄希破译细胞类型在其原生背景下的新状态。在这里我们向大家演示stLearn：一个综合分析以上三种数据类型的python库，stLearn首先像分析单细胞转录组一样识别细胞类型，不同的是stLearn可以在空间中重建组织内的细胞类型演变（Spatial trajectory inference ），并识别具有细胞间相互作用（Spatial cell-cell interaction）的区域。

stLearn的创新之处在于：

• SME normalisation ：(Spatial Morphological gene Expression normalisation，空间基因表达均一化)是一种组织内的均一化策略，它利用邻域信息(空间位置)和形态距离对基因表达数据进行均一化。试图把空间信息应用到数据分析的全流程。

• SMEclust：stLearn使用SME归一化数据，进行无监督聚类，将相似的点聚到聚类中，并根据组织中聚类的空间信息找到子类。与其他方法相比，SMEclust具有更高的空间聚类性能。

• spatial trajectory inference ：这部分是此演示文本的重点，我们将在下面描述之。

• spatial cell-cell interactions：为了研究细胞-细胞相互作用(CCI)，stLearn将细胞类型多样性和L-R共表达结合成一个统计量，可用于自动扫描整个组织切片。stLearn分别计算了L-R表达和细胞类型多样性，最后将这两种分析合并，以发现对细胞-细胞相互作用有重要作用的L-R对。

下面是stLearn分析流程的框架：

切入今天我们主要讨论的主题：空间轨迹推断

我们曾经表明：单细胞的一切分析 加前缀Spatial 都是一个新的分析点。trajectory 当然不会例外。

为了发现组织过程——例如，回答哪个癌细胞或克隆最先出现，或者癌症是如何进化的——stLearn提供了一种称为伪时空(pseudo-space-time, PST)轨迹分析的算法。PST是scRNAseq数据分析中常用的伪时间(pseudotime)概念的一个扩展，它被设计用来检测生物过程，这些生物过程可以从跨组织转录状态的梯度变化来推断。在PST中，进化轨迹是基于组织内细胞的空间环境和转录组图谱重建的。在SMEclust分析之后，我们实现了PST算法来寻找组织局部层次(即单个亚群的子群之间的关系)和全局层次(即亚群之间的关系)，分别对应亚群内部和亚群之间的关系。

stLearn首先使用基于全组织SME 均一化基因表达数据的PAGA轨迹分析，用于发现亚群内的联系。然后，利用稳健的轨迹推断方法——扩散伪时间(diffusion pseudotime ， DPT)方法计算伪时间。DPT方法使用类似随机游走的方法来测量细胞到细胞的转移。对于轨迹的方向，用户可以根据组织中正在研究的生物过程来定义根节点。然后结合基因表达值和物理距离计算伪时空距离( pseudo-space-time distance，PSTD)。在此基础上，构建有向图，像monocle一样寻找分支的方法类似：利用有向最小生成树算法对图进行优化，找出图中连接节点最短有根树和分支。可见，空间轨迹推断也是一种排序分析，只是构建距离矩阵的方法不同，这里的距离用到了空间信息。

工具安装：https://stlearn.readthedocs.io/en/latest/installation.html

数据下载：https://support.10xgenomics.com/spatial-gene-expression/datasets/1.0.0/V1_Adult_Mouse_Brain

安装分析工具和下载示例数据之后我们开始做演示：

import stlearn as st
from pathlib import Path
import scanpy as sc
import matplotlib.pyplot as plt
from matplotlib import cm as cm 
import numpy as  np
import matplotlib as mpl

mpl.rcParams["font.sans-serif"] = ["SimHei"]
mpl.rcParams["axes.unicode_minus"] = False
st.settings.set_figure_params(dpi=120)
# Reading data
data = st.Read10X(path="D:\\spatial\\V1_Adult_Mouse_Brain_filtered_feature_bc_matrix\\")

data
AnnData object with n_obs × n_vars = 2698 × 31053
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow'
    var: 'gene_ids', 'feature_types', 'genome'
    uns: 'spatial'
    obsm: 'spatial'

# 可以查看数据
data.obsm['spatial'] # 略
data.var['genome'] # 略 

可见 stlearn 用的是和scanpy一样的AnnData数据结构，其实stlearn的很多函数也是调用scanpy的，所以scanpy的所有分析在这里也是完全做的，如：

sc.pl.highest_expr_genes(data, n_top=20)

基本上与单细胞预处理一样的过程，质控，均一化，标准化过程。

# Save raw_count
data.layers["raw_count"] = data.X
# Preprocessing
st.pp.filter_genes(data,min_cells=3)

st.pp.normalize_total(data)
st.pp.log1p(data)

data.layers["normal_count"] = data.X
# Keep raw data
data.raw = data
data
st.pp.scale(data)
# 再存一份scale的数据
data.layers["scale_count"] = data.X

# 也存了三分数据：
print(data.layers['raw_count'])
print(data.layers["normal_count"])
print(data.layers['scale_count'])

降维聚类：

st.em.run_pca(data,n_comps=50,random_state=0)
st.pp.neighbors(data,n_neighbors=25,use_rep='X_pca',random_state=0)
st.tl.clustering.louvain(data,random_state=0)
sc.tl.umap(data)

data
AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

有了聚类信息可以看看分群情况：

labels = data.obs['louvain'].value_counts().index
students = data.obs['louvain'].value_counts()
colors =  ["#8da0cd","#fc8d62","#66c2a5","#ff7f00"]
explode = [0.1,0.1,0.1,0.1]

plt.pie(students,
        labels= labels,
        startangle=45,
        autopct="%3.1f%%",
        shadow= True,
        pctdistance= 0.7,
        labeldistance= 1.2 )
plt.title("cluster")
plt.show()

sc.pl.umap(data, color='louvain', palette=sc.pl.palettes.default_20)

st.pl.cluster_plot(data,use_label="louvain",tissue_alpha=1,spot_size=5,show_legend=True)

只画几个亚群：

st.pl.cluster_plot(data,use_label="louvain",tissue_alpha=1,spot_size=5,list_cluster=[1,2,3],show_legend=True)

绘制基因表达量：

st.pl.gene_plot(data,genes=["Cnp"],use_label="louvain",use_raw_count=True)

在特定亚群的表达情况：

st.pl.gene_plot(data,genes=["Cnp"],use_label="louvain",use_raw_count=True,list_clusters=[6,7])

基本的数据探索之后，我们开始做空间轨迹推断。注意，这里和单细胞转录组一样，根节点的选择也是很重要的。

import numpy as np
data.uns["iroot"] = np.flatnonzero(data.obs["louvain"]  == str(3))[50]
st.spatial.trajectory.pseudotime(data,eps=50,use_rep="X_pca")
data

AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain', 'sub_cluster_labels', 'dpt_pseudotime'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain', 'umap', 'louvain_colors', 'tmp_color', 'iroot', 'paga', 'louvain_sizes', 'diffmap_evals', 'threshold_spots', 'split_node', 'global_graph', 'centroid_dict'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca', 'X_umap', 'X_diffmap'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

没有空间信息的轨迹（拟时序）。

st.pl.non_spatial_plot(data,use_label="louvain")

st.pl.trajectory.pseudotime_plot(data,list_cluster="all",show_graph=True,node_alpha=1,tissue_alpha=1,edge_alpha=0.1,node_size=5)

st.spatial.trajectory.pseudotimespace_local函数可以构造局部轨迹, 该算法计算子类之间的时空距离。可以理解为亚群内部的轨迹（异质性）。

st.spatial.trajectory.pseudotimespace_local(data,use_label="louvain",cluster=3)
st.pl.subcluster_plot(data,use_label="louvain",cluster=3,tissue_alpha=1)

st.pl.trajectory.local_plot(data,use_cluster=3,branch_alpha=0.2,reverse=True)

可以看出指定的亚群内各部分间的转移轨迹，并给出对应的score值和方向。

欲构建全局的轨迹可以用：

st.spatial.trajectory.pseudotimespace_global(data,use_label="louvain",list_cluster=[0,3,4,5])

Screening:   1%|           [ time left: 00:15 ]Screening PTS global graph...
Screening: 100%|██████████ [ time left: 00:00 ]
Calculating: 100%|██████████ [ time left: 00:00 ]
Calculate the graph dissimilarity using Laplacian matrix...
The optimized weighting is: 0.47
Start to construct the trajectory: 3 -> 0 -> 5 -> 4

st.pl.cluster_plot(data,use_label="louvain",show_trajectory=True,list_cluster=[1,3,4,5],show_subcluster=False)

st.pl.trajectory.tree_plot(data)

有了在空间中的轨迹（形成不同的分支），我们肯定想知道哪些基因决定了这种轨迹，stLearn可以检测分支间过渡的marker gene。像差异表达分析的结果一样会返回一个gene list，这就可以走向基因调控，代谢通路之中了。

st.spatial.trajectory.detect_transition_markers_clades(data,clade=2,use_raw_count=True)
st.pl.trajectory.transition_markers_plot(data,top_genes=30,trajectory="clade_2")

data
AnnData object with n_obs × n_vars = 2698 × 18908
    obs: 'in_tissue', 'array_row', 'array_col', 'sum_counts', 'imagecol', 'imagerow', 'louvain', 'sub_cluster_labels', 'dpt_pseudotime'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'mean', 'std'
    uns: 'spatial', 'log1p', 'pca', 'neighbors', 'louvain', 'umap', 'louvain_colors', 'tmp_color', 'iroot', 'paga', 'louvain_sizes', 'diffmap_evals', 'threshold_spots', 'split_node', 'global_graph', 'centroid_dict', 'draw_graph', 'nonabs_dpt_distance_matrix', 'ST_distance_matrix', 'screening_result_local', 'PTS_graph', 'screening_result_global', 'clade_2'
    obsm: 'spatial', 'filtered_counts', 'normalized_total', 'X_pca', 'X_umap', 'X_diffmap', 'X_draw_graph_fa'
    varm: 'PCs'
    layers: 'raw_count', 'normal_count', 'scale_count'
    obsp: 'distances', 'connectivities'

data.uns['clade_2']

          gene     score       p-value
6842    Nap1l5  0.640180  1.287045e-60
6440   Tmem130  0.634474  2.987392e-59
2797     Rtl8a  0.632571  8.402284e-59
2659    Pcsk1n  0.630333  2.808432e-58
5676     Uchl1  0.629674  4.000150e-58
       ...       ...           ...
3480      Gmps  0.400526  3.161635e-21
9385    Slc1a6  0.400444  3.226337e-21
10329     Lsm6  0.400368  3.286873e-21
5349    Steap2  0.400231  3.399766e-21
13802     Gfap -0.546523  2.295220e-41

思考题：

考察下图，有了空间信息并检测出亚群边界后，这些边界附近的细胞有没有什么特殊的特征（处于边界本身就是一种特征），请尝试提出一些统计量来描述这些特征及其拟说明的生物学问题。

stLearn: integrating spatial location, tissue morphology and gene expression to find cell types, cell-cell interactions and spatial trajectories within undissociated tissues trajectories: Classes and Methods for TrajectoryData PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells Tempora: cell trajectory inference using time-series single-cell RNA sequencing data A comparison of single-cell trajectory inference methods Inference of multiple trajectories in single cell RNA-seq data from RNA velocity Untangling biological factors influencing trajectory inference from single cell data

---

References

[1] https://en.wikipedia.org/wiki/Trajectory_inference
[2] https://www.youtube.com/watch?v=Fbd08bIv4yk
[3] https://github.com/mdozmorov/scRNA-seq_notes
[4] https://stlearn.readthedocs.io/en/latest/Pseudo-space-time-tutorial.html

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注